May 1989 


LIDS-P-1874 


SIMULTANEOUS LINEARIZED INVERSION OF 
VELOCITY AND DENSITY PROFILES 
FOR MULTIDIMENSIONAL ACOUSTIC MEDIA 

by 

Ali Ozbek 

Department of Electrical Engineering and Computer Science 
and 

Earth Resources Laboratory 
Massachusetts Institute of Technology 
Cambridge, MA 02139 

Bernard C. Levy 

Department of Electrical Engineering and Computer Science 
University of California 
Davis, CA 95616 


This work was supported by the Air Force Office of Scientific Research under Grant No. 
AFOSR-85-0227, by the National Science Foundation under Grant No. ECS-87-00903, and 
by the Vertical Seismic Profiling Consortium at the MIT Earth Resources Laboratory. 



Report Documentation Page 

omZXZmss 

r • 1 1 r 1 11 • f r • • • ... . . 

VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for 

r any other aspect of this collection of information, 

1215 Jefferson Davis Highway, Suite 1204, Arlington 
failing to comply with a collection of information if it 

1. REPORT DATE 

may 1989 2. REPORT TYPE 

3. DATES COVERED 

00-05-1989 to 00-05-1989 

4. TITLE AND SUBTITLE 

Simultaneous Linearized Inversion of Velocity and Density Profiles for 
Multidimensional Acoustic Media 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

5c. PROGRAM ELEMENT NUMBER 

6. AUTHOR(S) 

5d. PROJECT NUMBER 

5e. TASK NUMBER 

5f. WORK UNIT NUMBER 

7. PEREORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

Massachusetts Institute of Technology,Laboratory for Information and 
Decision Systems,77 Massachusetts Avenue,Cambridge,MA,02139-4307 

8. PEREORMING ORGANIZATION 

REPORT NUMBER 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

10. SPONSOR/MONITOR’S ACRONYM(S) 

11. SPONSOR/MONITOR’S REPORT 
NUMBER(S) 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

Approved for public release; distribution unlimited 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

15. SUBJECT TERMS 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 

18. NUMBER 19a. NAME OE 

a. REPORT b. ABSTRACT c. THIS PAGE 

unclassified unclassified unclassified 

42 





Ozbek and Levy 


Velocity and Density Inversion 


2 


ABSTRACT 

The multidimensional inverse scattering problem for an acoustic medium is con¬ 
sidered within the homogeneous background Born approximation. The medium is 
probed by wide-band plane wave sources, and the scattered field is observed along 
straight-line receiver arrays. The objective is to reconstruct simultaneously the 
velocity and density profiles of the medium. The time traces observed at the re¬ 
ceivers ?ire appropriately filtered to obtain generalized projections of the velocity 
and density scattering potentials, which are related to the velocity and density per¬ 
turbations of the medium with respect to their nominal values. The generalized 
projections are weighted integrals of the scattering potentials; in two dimensions 
the weighting functions are concentrated along parabolas, in three dimensions they 
are concentrated over circular paraboloids. The reconstruction problem for the 
generalized projections is formulated in a way similar to the problem of x-ray, or 
straight-line tomography. The solution is expressed as a backprojection operation 
followed by a two or three dimensional space-invariant filtering operation. In the 
Fourier domain, the resulting image is a linear combination of the velocity and den¬ 
sity scattering potentials, where the coefficients depend on the angle of incidence of 
the probing wave. Therefore, two or more different angles of incidence are necessary 
to reconstruct the velocity and density scattering potentials separately. 
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INTRODUCTION 

In this paper, we consider the multidimensional inverse scattering problem for 
an acoustic medium within the homogeneous background Born approximation. An 
acoustic medium is probed by wide-band plane wave sources, and the scattered 
field is observed along straight-line receiver arrays. The objective is to reconstruct 
simultaneously the velocity and density profiles of the medium. 

The 1-D velocity/density profile inversion problem has been studied by a number 
of researchers, including Raz (1981a), Coen (1981), Hooshyar and Razavy (1983), 
Yagle and Levy (1984). This problem can in principle be solved exactly. 

The multidimensional problem, where the velocity and density profiles are al¬ 
lowed to vary in two or three dimensions, has also interested several researchers. 
For the multidimensional problem, no exact solution exists. Several approximate 
inversion techniques have been proposed, which linearize the scattering integral 
equations by using the Bom or Rytov approximations. In this context, it was 
shown that the experimental requirements of single parameter inversion problems, 
where the medium density is constant and only the velocity varies, and of mul¬ 
tiparameter inversion problems, where both the density and the velocity need to 
be reconstructed, are different. For the single parameter case, a single scattering 
experiment is sufficient to - at least partially - reconstruct the object of interest, 
whereas several experiments are necessary for multiparameter problems. 

Raz (1981b), and Clayton and Stolt (1981) have solved the multidimensional 
problem for an experimental setup where sources and receivers are available at all 
points on a plane, and scattered waves are measured at all frequencies. Coen, 
Cheney and Weglein (1984), and Ramm and Weglein (1984) have solved the same 
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problem for a complete set of sources and receivers, but using only two temporal 
frequencies. Hooshyar and Weglein (1986) have used two wide-band point sources, 
where the distance between the sources is required to be small compared to the 
distance from the sources to the scatterers. Beylkin and Burridge (1987) have 
presented a solution for a medium with a variable background, with sources and 
receivers surrounding the medium. 

The class of problems where the medium is probed by plane waves has been 
investigated by Norton and Devaney. Norton (1983) has used a flat transducer as a 
source of wide-band, plane wave illumination, and as a receiver of the backscattered 
waves. A second transducer, oriented at a different angle with respect to the first, 
is used as a receiver only. The two transducers are rotated together 180° around 
the object, and the scattered waves are recorded at all angles during the rotation. 
Devaney (1985) has extended the diffraction tomography theory to the variable 
density case. In this work the transmitted waves are measured on a plane parallel 
to the incident plane wave front, and the experiment is repeated for all view angles. 
In this scheme, two temporal frequencies are used in the insonifying wave. 

The present paper is a generalization of a previous work of the authors (Ozbek 
and Levy, 1987), where the multidimensional inverse scattering problem for a con¬ 
stant density acoustic medium was formulated and solved as a generalized tomo¬ 
graphic problem. In this paper, we similarly filter the time traces observed at the 
receivers to obtain generalized projections of the velocity and density scattering po¬ 
tentials, which are related to the velocity and density variations in the medium. The 
generalized projections are weighted integrals of the scattering functions; in the two- 
dimensional geometry the weighting functions are concentrated along parabolas, in 
the three-dimensional geometry they are concentrated over circular paraboloids. 
Thus the inverse scattering problem is again posed as a generalized tomographic or 
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integral geometric problem. 

The reconstruction problem for the generalized projections is formulated in a 
way similar to the problem of x-ray, or straight-line tomography. The solution 
is expressed as a backprojection operation where we sum the contributions of all 
projections passing through a given point in space, followed by a two or three 
dimensional space-invariant filtering operation. In the Fourier domain, the resulting 
image is a linear combination of the velocity and density scattering potentials, where 
the coefficients depend on the angle of incidence of the probing wave. Therefore, 
two or more different angles of incidence are necessary to solve for the velocity and 
density scattering potentials separately. 

The main difference between the approach that we propose here and the diffrac¬ 
tion tomography technique developed by Devaney (1985) is that we use wide-hand 
plane waves at just a few angles of incidence to reconstruct the medium density 
and velocity, whereas the diffraction tomography formulation relies on narrowband 
plane wave sources at just two, and in practice several, frequencies, but for all 
angles of incidence. These two approaches are in some sense dual of each other, 
since they trade wavevector orientation against frequency. This distinction leads to 
significantly different processing algorithms, and in fact, as mentioned above, the 
inversion procedure proposed here is significantly closer to x-ray tomography than 
to diffraction tomography. 

Also, the different choice of experimental conditions appearing here is a reflection 
of the difference existing between medical imaging applications, which motivated 
the diffraction tomography approach, and geophysical tomography problems, which 
are at the origin of the present work. In medical imaging it is possible to rotate the 
object to be imaged, i.e. the patient, over a 360° range, but in exploration geophysics 
only a very limited range of angles of incidence can be achieved for surface, or even 
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vertical seismic profiling recordings. 

The paper is organized as follows: we treat the 2-D case in detail in Sections 
I-III, and just summarize the results for the 3-D case in Section IV. In Section I, 
the inverse scattering problem is formulated within the Born approximation and 
redefined as a generalized tomographic problem. The backprojection operation is 
defined and related to the generalized projections in Section II. In Section III, the 
separate reconstruction of the velocity and density scattering potentials is discussed. 
In this context, it is shown that substantially more than two plane-wave experiments 
are required in order to be able to recover the velocity and density perturbations 
in a numerically reliable way. We summarize the results for the 3-D geometry in 
Section IV. A 2-D numerical example is presented in Section V, and Section VI 
contains some conclusions. 

L PROBLEM DESCRIPTION 

In this paper we will treat the two-dimensional case in detail, and summarize the 
results for the three-dimensional case. Consider the scattering experiment described 
in Fig. 1. A 2-D acoustic medium is probed by a wide-band plane wave and the 
scattered field is observed along a straight-line receiver array. The Fourier transform 
P(x, w) of the pressure field at point x= {x,y) satisfies (Chernov, 1960) 

‘ 

where c(^ is the propagation velocity, and p(z) is the density of the medium at 
point X. The acoustic equation (1) can be rewritten as 


(V^ + k^)P{x,uj) = -k^U,{x)P{x,uj) -F • VP(x,w), 


(2) 
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where 


UM) ^ 


c^x) 


U^^x) 4 In^, 
Po 


(3) 

(4) 


and k = ujco is the wavenumber. Here, £7c(x) and are respectively called 

the velocity and the density scattering potentials. Co and Po are respectively the 
velocity and the density of the background medium. We assume that c(i) and p(^ 
do not deviate significantly from their nominal values c<, and po\ consequently 14 (x) 
and Up{^ are small with respect to 1. We also assume that 17c(x) and f4(x) have a 
bounded support V , which is located completely on one side of the receiver array. 

The probing wave Po{x,u}) satisfies 


(V2 + A:2)P,(^a;) =0, 


(5) 


so that the scattered field P,(x, w) = P(x, w) — Po(^ u) obeys 

(V^ + fc2)P,(x,a;) = -k^Uc{x)P{x,uj) + VU,{x) • VP{x,u). (6) 

The solution of (6) is given by the Lippmann-Schwinger equation (Taylor, 1972) 

Ps{x, uj)= j dx' [A:'l7e(x')P(x', w) - V17p(x') • VP(x', cj)] Go{^ x', u), (7) 

where G(,(x,x',a;) is the Green’s function eissociated with a point source in a homo¬ 
geneous medium: 

(V* + k^)Go{x,x!,u}) = -d(x - x'). (8) 

Equation (7) demonstrates the nonlinear relation that exists between the po¬ 
tentials £4(x), £4(x) and the pressure field P(^w). To linearize this equation we 
adopt the Born approximation, whereby we assume Ps{x,u) C Po{^u}). Then the 
Lippmann-Schwinger equation becomes 

Ps{x,uj) « I dx' [k^Uc{x;)Po{x!,u) - VC^,(x') . VP<,(x',a;)] G<,(x,x',a;). (9) 
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The Born approximation assumes that the scattered field P, (x, w) is small through¬ 
out the volume of the object, which requires that both the magnitude of the scat¬ 
tering potentials Ue and Up be small and that the volume of the region V be small 
with respect to the dominant wavelength of the probing wave. 

If, instead of the Born approximation one uses the Rytov approximation, the 
requirement that the size of the scattering region be small can be relaxed (Chernov, 
1960; Tatarski, 1961; Devaney, 1981). The Rytov approximation is obtained by 
representing the total pressure field P{^u) in terms of its complex phase and lin¬ 
earizing the resulting Riccati equation satisfied by the phase fiuctuation (Devaney, 
1985). 

On the other hand the Bom approximation is more accurate than the Rytov 
approximation for refiected waves (Beydoun and Tarantola, 1986). For the setup 
considered in this paper, where the wideband property of the probing wave replaces 
the large number of view angles available in diffraction tomography, it was shown 
in Ozbek and Levy (1987) that the reflected wave configuration provides the most 
coverage in the Fourier domain for a bandlimited source. Therefore, we adopt the 
Born approximation in this paper, although similar results can also be derived for 
the Rytov approximation. 

Next, we simplify the second term in the integrand in equation (9) by using the 
identity (Norton, 1987) 

{VUp • VP,)Go = V ■ (UpG.VPo) - UpV • [GoVPo), (10) 

and applying the divergence theorem. Over a surface S located outside the domain 
T where the density inhomogeneities are located, we have Up{^ — \n{p{^/po) = 0, 
so that 

dx'V . [Up{P)G,{x,^,u})VP,{^,uj)] 



Ozbek and Levy 


Velocity and Density Inversion 


= j dsUp ( 2 ')Go(x, s', w)n{ 5 ) • VPq( s', w) 

= 0 . ( 11 ) 

For the second term in equation (10), we have 

V • (GoVPo) = VGo • VPo + GoV^Po = VGo • VPo - k^GoPo- (12) 
The incident wave is given as 

PoiP,u) = e'’‘^-^', (13) 

where £ = (cos0, sin 5) is the unit vector which indicates the angle of incidence 
of the plane-wave source. Therefore, within the Born approximation, the scattered 
field at a receiver point ^ is given by 

+?7,(x')V,.Po(x',u;) . V..Go(i,x',u;)} 

= J dx' |A:^[17 c(x') — ?7p(x')]Go(|^, x', oj) 

+tA:l7,(x')[£- V,.God,x',w)]} Po(x',a;). (14) 

To compare the far-field scattering patterns that are due to velocity and density 
perturbations, let us use the far field approximation k\^ — £j » 1. We find 

V*»Go(£,x',a;) « —||■God,x',u;), (15) 

12 . - £1 

which is valid for both the 2-D and 3-D Green’s functions. Then (14) becomes 

Po(e,w) « k^l d^{u,{^) - [1 + cosad,x',£)]C/,(x')} Go(e,£', w)Po(x',a;), (16) 

where a(£, x', £) is the angle between the vectors —£ and ^ — x'. From a physical 
point of view, a:(^,x',£) is the angle at point x’ between the ray linking x' to the 
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soiirce, which is specified by vector —£ since the incident wave is a plane wave, 
and the ray linking z' to receiver Since the background medium is assumed to 
be constant, these rays are straight lines. The different weights in front of i7c(x) 
and f7^(z) in (16) correspond to different scattering patterns which are plotted in 
Figs. 2a and 2b. As indicated by Fig. 2a, the scattering pattern due to is 

that of a monopole, whereas Fig. 2b shows that the scattering pattern due to Up{^ 
is that of the sum of a monopole and a dipole. Therefore, the scattering due to 
density perturbations is most prominent for reflected waves, and least prominent 
for transmitted ones. 

Equation (14) can be written in a slightly different form if we introduce the 
scattering potential 

(17) 

associated to the compressibility /c(z^) = l/c^(^/)(x). We find 

( 18 ) 

for c(x) near Cg, as assumed above. Therefore (14) can be written as 

P,{(,u) = J V,.G„(f,i,',w)]} P<,(£',w), 

(19) 

and the objective of this paper will be to solve this integral equation for the com¬ 
pressibility and the density scattering potentials C4(^ and ?7p(x), or equivalently 
for f/p(x) and the velocity scattering potential ?7c(x). 

For the 2-D geometry under consideration, the Green’s function is given by 

<^o(£>£.*>w) = ^H^^^(fc|x — x'l), (20) 

where H[)^^(') indicates the Hankel function of order zero and type one. In the 
following, it will be assumed that the receivers are located along a straight line 
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perpendicular to the unit vector 0 = (cos^, sin^) and whose distance from the 
origin in the direction ^ is p, as shown in Fig. 1. The position of an arbitrary receiver 
along this line is therefore given by ^ = p^ + , where 0'*’ = (sin — cos is a 

unit vector perpendicular to and ^ is an arbitrary coordinate. Then (19) can be 
expressed as 

( 21 ) 

where 

F.U.k) = '-^1 ixV.(£')e“2^'HS,‘>(J:|s' - el), (22) 

FM, *’) = ^ / <ix'Cl,(x')e“i'^' {« • V.,H<‘>(l=|s' -11)} 

= -|/<fa'P.(£')«“2-'[£-j|^]H<‘>(fc|x'-e|), (23) 

and where = —dIl^Q\v)/dv is the Hankel function of order 1 and type 1. 

We define the inverse Fourier transform of F{^,k) with respect to k as p(^, r): 

SU, r) ^ i dkp{t 4)«-“^ (24) 


»■) and pp(^, r) are similarly defined as the inverse Fourier transforms of .F«(^, k) 
and Fp{^,k), respectively. Taking into account the fact that (see Morse and Fesh- 
bach, 1953, pp. 1362-1363) 


7''{fHW(4a)} 

l(r-u) 

Vr2-u2’ 

(26) 

where !(•) is the unit step function, and the relationship 


{-|H<‘>(iu)} 

_ r l(r - u) 
u y/r^ - u2 

(26) 

that follows from (25), we find that 




1 -ffp(^,r). 

(27) 
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where 


and 


S.(e.r) = / d^u,[x’) 


l(r - £ • x' - |x' ^|)_ 




= j d^U,{^) S 


- i)l I' - i- £'l H' -i-s' -Is' -II) 

Is' - ii J [ w - ii J v'(>--£•$')*-Is' -ir 


(28) 


(29) 


Equation (28) expresses as a weighted integral of the compressibility 

scattering potential ^/^(x) where the weighting function is nonzero in a region with 
parabolic support. The parabola satisfies the equation r = £• x — |x — £| where r, 
£ and £ are given and x varies. The directrix of this parabola is the line £ • x = r 
which is perpendicular to the direction £ of incidence of the probing wave, and 
whose focus is the receiver point The weighting function becomes infinite for 
values of x along this parabola, so that the largest contribution to the integral is 
made by the values of 17,t(x) which lie along the parabola. In some sense, SfK(£,r) 
is a projection of U^ig^ with respect to a function whose singularities are algebraic 
and located along a parabola. 

In (29), gp{^,r) is similarly expressed as a weighted integral of the density scat¬ 
tering potential The weighting function again has a parabolic support, but 

it contains two additional factors. The first factor, £ • (x' — £)/Ix' —^|, is again 
identified as cosa(^,x',£), where o:(£,x',£) is the angle between the vectors —£ and 
£ — x'. The second factor in (29), (r — 0- x!)/\x' — £|, is a term which equals unity 
on the parabolic wavefront of the weighting function, and grows as the focus of the 
parabola (the receiver point) is approached. 

In the following, it will be assumed that the projections g{$,r) constitute the 
data obtained from a single plane-wave scattering experiment. To see why this is 
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the case, observe from (21) and (24) that 

Then since 

/««)} = ^ + (31) 

and observing from (21) that P,(£,a; = 0) = 0, we find that 

= -27rCo I / ^ dT f dsP,{^,s) - ^ f dr f dsP,(|_,5)|. (32) 

00 j—OO Z j--00 — oo ) 

The projections sf(£, r) are therefore obtained by integrating twice the time domain 
scattered field Pj(^, i) observed at and then subtracting a constant equal to half 
the value of the double integral at i = +oo. This shows that the knowledge of the 
observed scattered field P,(^, t) for all t is equivalent to that of projections g{$, r) 
for all r. 

On the basis of the above observations, the inverse scattering problem can be 
formulated as follows: given the generalized projections 

{&(^)'’) : —oo < ^ < oo, —OO < r < oo}, 

we want to reconstruct the scattering potentials Uk,{^ and 17^(x). 


II. BACKPROJECTION OPERATION 


Proceeding as in the constant density inversion problem treated in Ozbek and 
Levy (1987), the first step of our inversion procedure is to perform a backprojection 
operation on the projections Sf(^,r). We define it as 


Ub{x) = 


Ijr-i-x- |x-||) 
v/(r-£'z)2 _ 


(33) 
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At a given point x, this operation sums the contributions of x to the projections 
g{^,r), as appearing in equations (27)-{29). By performing this backprojection 
operation for every point in the plane, this gives an image, 17 b (x). 17 b(x) is the 
analog of the backprojected image obtained in x-ray tomography (Deans, 1983) by 
summing the line projections passing through a given point. 

Our first objective is to relate Ub{^ to the projections g{^,r) in the frequency 
domain. It can be shown that the 2-D Fourier transform of 17 b (x) is given by (see 
Ozbek and Levy, 1987, Section 3 and Appendix A) 

UbOc) = I dxUBi^e-'^-^ 


tit 

k7l 





A • ■0^ 
2&J ’ 


K 



(34) 


where 


Uk(, kr) = r r drg{$, (35) 

J-OO J-OO 

is the 2-D Fourier transform of sr(^,r), k= {kx,ky), k = |fc|, A = (fc^ - ik®, 2kj,ky), 
0 = (cos(0 -I- 0), sin((9 + 0)), and = (sin(0 -1- 0), - cos(0 -|- 0)). 


III. SEPARATE RECONSTRUCTION OF ft(i) AND Up{k) 

In this section, we first derive a frequency domain relationship between the pro¬ 
jections g(^,r) and scattering potentials C4(x) and 17p(x), thus obtaining a ’’Pro¬ 
jection Slice Theorem” (Deans, 1983) associated with the problem. Inverting this 
relationship provides both a frequency domain relationship between ^,)(i), 

and Ub{}^, and a method for the separate reconstruction of 17„(x) and 17p(^, or 
equivalently, of 17«(^ and 17p(x). 

From (27), we have 


Uk(,kr) = h{k(,kr) - 'g^{k(,kr). 


( 36 ) 
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where kr) and g^(k(, kr) are the 2-D Fourier transforms of O and gp{^, r), 

respectively. Since ^^(^ 5 r) = g{^,r) for the constant density case, it was shown in 
(Ozbek and Levy, 1987, Section 4 and Appendix B) that gi^{k^,kr) can be expressed 
as 


L(*s. kr) = cr, = krl+ + ^ffi(kr)slkl - k\ ^ 


(37) 


for \k(\ < |A:r|, where 

UM = / dx't7,(x')e-‘^' (38) 

is the 2 -D Fourier transform of Uk.{^, and 

^ j - 1-1 if X • p > 0 for all X G "V, 

(-1 if X • p < 0 for all X € V. 


7 describes which side of the receiver array the support V of the inhomogeneities is 
situated. 

We now express p^(/:f,A:,.) as a function of the 2-D Fourier transform Up{h) of 
f7p(x). We first take the Fourier transform of gp{^^r) with respect to r: 


Mi<kr) = j^^drg,[(,r)e 
= F;{(,kr) 

" Wr I {a- -Q)}, (39) 

where Fp denotes the complex conjugate of Fp and Ho*^(-) is the Hankel function 
of order 0 and type 2. Now, taking the Fourier transform with respect to ^ gives 


Uk„kr) = i/ d^Url,^)e-<-'>---{i-V.jN{x:,k„k.)}, 


(40) 
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where (Ozbek and Levy, 1987, Appendix B) 

N{^,k^,kr) = 

J-oo ■" 

= 2sgn(<:,) 

for Ifcfl < \kr\. Consequently, we get 

•Uf,(jc- krt + + 7Sgn(A:r)^fcJ - ^ , (42) 

for \k^\ < jfcrl. Combining (36), (37) and (42), we obtain 


kr) 


- 5 ^ e<- {U. (i= a+ ..r + 7»gn(*,)VSm) 

i- [k^$^ + 7sgn(fc,)y^A:2 - . .x __\1 

-- Uf, (^k = kri. + k($ + ')sgn{kr)yjk^ - k^ i, 


(43) 


for |A:f| < |A:,.|. For |A:j| > |fc,.|, g(^k^,kr) is related to the part of the observed 
scattered field that corresponds to evanescent waves (Ozbek and Levy, 1987), and 
we do not make use of this portion of g{k^, kr) in our inversion scheme . The inverse 
formula of (43) is 


Ur{!^ = 


UM - 2 U,{k) 

27rk -e ^ ^ 



2k-i 



keC, (44) 
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where the cone C is defined below, and where Uc{k) = UkIM} + Up{^ is the 2 -D 
Fourier transform of the velocity scattering potential Ue{^. Here, denotes the 

2-D Fourier transform of the reconstructed potential Uji{^ obtained by applying 
the constant density reconstruction procedure to the projections ff(^,r) obtained 
from a variable density and velocity medium. 

Equation (43) represents the "Projection Slice Theorem” associated with the 
variable density inverse acoustic problem relating the 1-D Fourier transform of 
^r) with respect to ^ to a semicircular slice of the 2-D Fourier transform of 
17 r(^. For a fixed k^, ff(k(,kr) gives Ur(!^ along a semicircle of radius |Ar,.| centered 
at as shown in Fig. 3. By letting kr vary, these semicircles span a cone C, which 
is defined as 

C = {k: \k-P\>V2/2} (45) 

for 7 = +1, where {k,k) and (cos((5 + ^)/2), sin((^-|-0)/2)). The angular 
range of this cone is 90°, as indicated in Fig. 3. For 7 = — 1 , C is the complement 
C of the above cone. 

Combining (34) and (44) gives 

UR{k) = ^^^UB{k) = U4k)-2cos\U,{k), kec. (46) 

where ^ is the angle between the vectors k and £, with k= (k,k). This relation is 
the key result of our paper. It shows that the reconstructed image UR(k) can be 
obtained by applying the 2-D filter 7 A • ^/27r^ to the backprojected image t/slk). 
This relation is similar to the identity which is used for 2-D backprojection and 
filtering x-ray reconstruction methods (Deans, 1983). Since the filter A • ^ diverges 
for high frequencies k, it needs to be "clipped” as k becomes large. Furthermore, 
identity (46) shows that Ur(}£^ is a linear combination of Uc(J^ and Up{k). This 
implies that it is not possible to reconstruct these two potentials from a single 
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experiment. To reconstruct Ue{k) and Up{k) separately, we need in principle two 
experiments with, plane waves incident at angles 0^ and 0^] then we can solve the 
system 


1 -2{k'L? 

um 


UrM 

1 

. UAH) , 


UR2ik) 


which requires inverting the matrix 

For the numerical stability and robustness of the matrix inversion procedure, 
the matrix M{^0i,0^) must be as nonsingular as possible. The most appropriate 
measure of the singularity of a matrix is the smallest singular value of the matrix. 
The smallest singular value of the matrix M is 

= nhna,- such that det(o-?J — M'M) = 0 

= 1 - [(277^ + 2r)t + l)' - 4 {n\ - , (48) 

where tji = fc • and 772 = 

Inversion of M would be most robust when OxmvXM) maximized. This takes 
place for values of and 0^ such that = 0 and % = ±0i or = ±0^,. 

Therefore the two probing waves must be incident at angles perpendicular to each 
other. Under this condition, let us consider the frequency domain coverage we 
would have for finite bandwidth data, assuming that we have receiver coverage 
surrounding the medium. Neglecting the low frequency cutoff band, the frequency 
domain coverage due to a single probing wave has a "figure-of-eight” shape aligned 
with the direction of the probing wave (Ozbek and Levy, 1987). When two probing 
waves are used, I7(.(^) and Up{!^ can be solved only in regions where where there is 
double coverage, as indicated by the shaded areas in Fig. 4. However, if we consider 
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the superimposed "radiation pattern” of <Tniin(M) drawn in Fig. 4 also, we see that 
M is most singular for those values of k where we have double coverage. 

For general values of and 0^ the situation is similar. In general, M is singular 
for values of which satisfy • £i| = |^ • ^|. Therefore, for ^ = ±^, M is singular 
forall^ otherwise, it is singular for fc = ±(ii+l 2 )/l^+^l orfc = ±(ii—£ 2 )/lii“i 2 l- 
These are the directions which in fact bisect the regions where there is double 
coverage. The coverage and radiation pattern for the case when the angle between 
and 0^ is 45° is shown in Fig. 5. 

The preceding analysis confirms what researchers in diffraction tomography and 
exploration geophysics have suspected for some time: namely, that it is exceedingly 
difficult to reconstruct simultaneously the velocity and density perturbations of 
an acoxistic medium in a numerically stable way. A purely theoretical analysis 
neglecting the bandlimitation of the probing wave would lead to the conclusion that 
only two experiments at different angles of incidence are necessary. Yet, as shown 
above, the numerical results obtained by such an approach would be worthless. In 
practice, one must use substantially more than two angles of incidence, say angles 
0J., • • •»£iv> for each k, solve the resulting system 


1 -2(k-ij^ 

1 -2(k-lj^ 


' Uc{k) 



.1 -2{k-ij\ 


. U,{k) 




M(k) 


by the least squares method, where {ty, ..., 4i>} C {1,2, ...,iV} is the set of 
indices corresponding to the angles of incidence for which the probing wave provides 
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coverage at k. This gives the solution 

UM 


(50) 


IV. THREE-DIMENSIONAL GEOMETRY 

After discussing the 2-D experimental geometry, we summarize the correspond¬ 
ing results for the 3-D case. For the 3-D geometry, we assume that the receivers are 
on a plane; for convenience we choose this to be the x-y plane. The position of an 
arbitrary receiver located in this plane is therefore given by £ = (^, 0), where ^ 
represents the x-y coordinates of the receiver. The Green’s fimction due to a point 
source is 


Co{x,g^,u)) = - 




(51) 


A'n\x — ^\’ 

Under the Born approximation, the Lippmann-Schwinger equation takes the form 
(19), and the projections in this case become 


= &«(|j„r)-y^(|y,r). 


where 


(52) 


(53) 


and 


k'-ii 


%'-ir 


( 54 ) 
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where £ is the unit vector in the direction of propagation of the incident wave. The 
projection r) is again a surface integral of the scattering potential over a circu¬ 

lar paraboloid, like in the constant density case. The weighting function appearing 
in the representation (53) of is an impulse, and in this sense the 3-D case 

is quite different from the 2-D case. Projection gp{$^,r) shows further deviations. 
Like in the 2-D case, it has a cosine dependence on the angle between the direction 
of the incident wave and the direction of scattering. The weighting function, in 
addition to an impulse, contains a smoother term which in fact makes gp{^^,r) a 
noncausal function of r. This is due to the l/k^ filter that is applied to the scat¬ 
tered field P,(£, cj) to obtain the projection grp(|^,r). If the far-field/high-frequency 
approximation A:|x' — 1 is made, the impulsive term clearly dominates. 

We introduce the backprojection operation as in the constant density case: 

Us{^ ^ jdi ^/_” (55) 

In the frequency domain, the backprojected image can be related to the parabolical 
projections as (see Ozbek and Levy, 1987, Section 8 and Appendix C) 



where 

Ub{}^ = / dzUei^e-^^-^ 

gik^T^K) = y 

are the 3-D Fourier transforms of Ub{^ and ff(^y,r), and 

A, = {kl-kl- kl, 2k,k„ 2k,k,), 
A, = {2k,k„ kl~ kl-kl, 2k,k,). 


(57) 

(58) 
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From (52), we again have 


fCr) = L(k^T> M - ^r), (59) 

where and gp(k^Tykr) are the 3-D Fourier transforms of fl'K(£j., r) and 

ffp(ir^r), respectively. Again, since S'«(|y,r) = g(i^,r) for the constant density 
case, it was shown in (Ozbek and Levy, 1987, Section 8 and Appendix D) that 
ff)c(^rj^r) takes the form 

L(^r, kr) = - = k,i + (^r,7sgn(fc,)>/fcJ^^]fc^)) , (60) 

for l^rl < lA:^!, where U^nik) is the 3-D Fourier transform of 17^(^, and 

^ I -hi if 2 > 0 for all s 6 V, 

[ —1 if z < 0 for all X e V. 

Using the intermediate results derived in (Ozbek and Levy, 1987, Appendix D), 
we also have 


dp(k^T,kr) 


[^- (k4T,7sgn(kr)y^k^ - I^tH] 


Ikrly/k^ - I^tI^ 

•Uf [k = kri + (^r, isgn{kr)yjkl - |^rP)) , (61) 

for l^rl ^ lA^rl? where U^(^ is the 3-D Fourier transform of Up{^. 

Combining (52), (59) and (60) yields 

Uk^T, kr) = 

~ (^r,7sgn(M)//:2 - |^r|2)) 

£• ik^T,lsgn{kr)Jk^ - . / .. , _ 

- Up [k = krO -h (^r, ^sgn{kr)\/k^ - |^tP))} , 

(62) 
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for |^y| < \kr\. 

The inverse formula of (61) is 


c^*(a = 


i2k4sL-MM] 

47rk‘i [2k-d^2k-i] 



keC, 


(63) 


where 

= [2k^k,, 2kyk„ kl-kl- kl), 

and where UeiM) = ^(c(^) + ^#>(^) is the 3-D Fourier transform of the velocity- 
scattering potential l!7c(x). Here, t/n(fc) denotes the 3-D Fourier transform of the 
reconstructed potential Ur{o^ which is obtained by applying the constant density 
reconstruction procedure to the projections g{^, r) obtained from a variable density 
and velocity medium. 

For l^rl ^ l^r|> S'® in the 2-D case, 5(^r>^r) is related to the part of the 
observed scattered field that corresponds to evanescent waves, and we do not make 
use of this portion of g{k^Ti ^r) in our inversion scheme. 

Equation (61) represents the "Projection Slice Theorem” for the 3-D case as¬ 
sociated with the variable density inverse acoustic problem. For a fixed kr, the 
2-D Fourier transform of g(£y,fcr) gives the 3-D Fourier transform of Ur{^ over a 
hemisphere of radius |fcr| centered at krO. By letting kr vary, U{k) is determined in 
a cone C, which again covers half of the 3-D frequency space. 

Combining (56) and (62) gives 


Mi) = -2cos’(U,(k), keC. 


(64) 
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Therefore, Ur{!^ for fc € C can be obtained from Usik) by a 3-D space invariant 
filtering operation. 

UR{k) is related to Udb) and Up(k) with the same relationship ((46) and (64)) 
in two and three dimensions. Therefore, the procedure for individual reconstruction 
of Uc{^ and £/p(^ is the same in these two cases. 

V. NUMERICAL EXAMPLE 

The theory presented in this paper was tested for the two-dimensional case, using 
computer-generated synthetic data. Figs. 6 and 7 show the velocity and density 
scattering potential models, and 17p(x), respectively. The scattering potentials 
correspond to velocity and density anomalies which are constant in square-shaped 
areas of dimensions 35 m x 35 m. The background medium was homogeneous 
with a velocity of 5000 m/s and a density of 2000 kg/m^. The magnitude of the 
velocity and density perturbations were 10% of the nominal values. A lowpass 
source wavelet with a cutoff frequency of 425 Hz was used, so that the object sizes 
are three times the shortest wavelength in the source signal. The regions of anomaly 
are separated by a distance equal to six times the shortest wavelength. The synthetic 
scattered waves were obtained by using the forward scattering equation under the 
Born approximation; however, since the velocity and density perturbations are of 
limited size with respect to the shortest wavelength, and the magnitudes of these 
perturbations are relatively small, we do not feel that the approximation is critical 
for this particular example. The entire image area was 500 m x 500 m, the grid 
size was 5 m x 5 m, and receivers were located on all sides around the medium, 
100 on each side. 

As indicated in Section III of this paper, in order to guarantee the numerical 
stability of the procedure for reconstructing separately the velocity and density’ 
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inhomogeneities, more than two sources are needed. In this experiment, we have 
used eight angles of incidence, at 22.5® intervals. The inversion was performed over 
the regions in the k domain where coverage was provided by at least five probing 
waves; i.e., using the notation of eq. (49), N — S, P > 5, and rank(M) > 4 for all 
inversion points k. This corresponds to carrying out the inversion over a circular 
lowpass region with a radius of about 55% of the maximum frequency coverage 
provided by a single source. 

For comparison, we first examine images ?7 b"(z} and U^{^. To obtain these 
images, in the frequency domain, values obtained due to different sources providing 
multiple coverage were simply averaged point by point. Fig. 8 shows the back- 
projected image can be interpreted as a "migrated” image of the 

velocity field for a constant density medium (for migration, see, for example, Claer- 
bout, 1985). Fig. 9 depicts Z7^"(x), which is the image obtained by applying the 
constant density reconstruction procedure to the data obtained from a variable 
density and velocity medium. 

Some observations can be made regarding these images. Both images display the 
locations of the scatterers; however the "inversion” image is much better focused 
than the "migration” image. This effect has also been noted by other researchers 
(Esmersoy and Miller, 1987). In addition, the values of differ by orders of 

magnitude from the numerical values of the true scattering potentials. On the other 
hand, looks like Uc(^ — Up(x), and actual constructed values confirm this. To 

interpret this result, observe from (46) that, for the "ideal case" where plane wave 
experiments are performed for all angles 0 of incidence, the averaging scheme de¬ 
scribed above for combining the reconstructed images obtained for different probing 
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angles can be written as 

® = Ii' -Upi!s) = VM, (65) 

where Ur corresponds to the inversion result for one source, and corresponds 
to the result after angular averaging. Therefore, averaging the reconstructed poten¬ 
tial Ur{J^ over different angles is equivalent to reconstructing the compressibility- 
potential of the medium. 

Figs. 10 and 11 show the separate reconstructions of t/c(x) and respec¬ 

tively. The numerical values obtained are within 20% of the model values. There 
exist several sources of error. The first of these is the bilinear interpolation proce¬ 
dure which is used to convert the discretized 2-D Fourier transform k^) of the 
projections into the discretized Fourier transform Ur{!c) of the reconstructed image. 
A second source of error is the fact that all the inversion results developed in this 
paper assume that receiver arrays are infinite, whereas the arrays which are used 
for the present example have a finite length. Other errors are due to the fact that 
the source wavelet is bandlimited, and needs to be deconvolved. 

Finally, the DC levels of the velocity and density scattering potentials cannot be 
reconstructed separately with the derived inversion formulas, since the coefficient 
of Up{k) in equation (44) is not analytic around ^ = 0. The DC level of U^{x) can 
be recovered from equation (22), so that if the DC level of either the density or the 
velocity is known, the other one can be computed. In our implementation, we have 
estimated Uc{k. = 0) and Up[k = 0i) as a weighted average of the closest eight values 
in the discrete wavenumber domain. Adopting the reciprocal of the square of the 
distance as the measure of weight, we assigned a weight of 1/6 to the closest four 
samples, and 1/12 to the next closest samples which are diagonally located. 
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VI. CONCLUSIONS 

In this paper we have considered the problem of the reconstructing separately 
of the velocity and density inhomogeneities of a multidimensional acoustic medium 
probed by wide-band plane waves. The problem was posed as a generalized tomo¬ 
graphic problem, where weighted integrals of the velocity and density scattering 
potentials Uc(x) and Up[x) are used as data. A backprojection operator C/s{x) was 
defined, which was related to the generalized projections in the Fourier transform 
domain. It was shown that, by applying a time-invariant filter to 17 b (^, we can 
obtain an image, Ur{c^, which in the Fourier domain is a linear combination of 
the velocity and density scattering potentials, and where the coefficients depend 
on the angle of incidence of the probing wave. Therefore, for numerical stability, 
several angles of incidence were used to solve for the velocity and density scattering 
potentials separately. 


ACKNOWLEDGEMENTS 

This work was supported by the Air Force Office of Scientific Research under 
Grant No. AFOSR-85-0227, by the National Science Foundation under Grant No. 
ECS-87-00903, and by the Vertical Seismic Profiling Consortium at the MIT Earth 
Resources Laboratory. 



Ozbek and Levy 


Velocity and Density Inversion 


28 


REFERENCES 

G. Beylkin and R. Burridge (1987). ”Multiparameter Inversion for Acoustic and 
Elastic Media,” presented at the 57th Ann. Int. Mtg. and Expos., Soc. Explor. 
Geophys., New Orleans, Expanded Abstracts, 747-749. 

W. B. Beydoun and A. Tarantola (1986). "First Born and Rytov Approximations: 
Modeling and Inversion Conditions in a Canonical Example,” submitted to J. 
Acoust. Soc. Am. 

L. A. Chernov (1960). Wave Propagation in a Random Medium (Dover, New 
York). 

J. F. Claerbout (1985). Imaging the Earth’s Interior (Blackwell, Oxford). 

R. W. Clayton and R. H. Stolt (1981). ”A Born-WKBJ Inversion Method for 

Acoustic Reflection Data,” Geophysics 46, 1559-1567. 

S. Coen (1981). "Density and Compressibility Profiles of a Layered Acoustic 

Medium from Precritical Incidence Data,” Geophysics 46, 1244-1246. 

S. Coen, M. Cheney, A. Weglein (1984). "Velocity and Density of a Two-Dimensional 
Acoustic Medium from Point source Surface Data,” J. Math. Phys. 25, 1857- 
1861. 

S. R. Deans (1983). The Radon Transform and Some of its Applications (Wiley 
Interscience, New York). 

A. J. Devaney (1981). "Inverse Scattering Theory within the Rytov Approxima¬ 
tion,” Opt. Lett. 6, 374-376. 



Ozbek and Levy 


Velocity and Density Inversion 


29 


A. J. Devaney (1985). ’’Variable Density Acoustic Tomography,” J. Acoust. Soc. 
Am. 78, 374-376. 

C. Esmersoy and D. Miller (1987). "Stacking Versus Back Propagation in Seismic 
Imaging: Duality for Multidimensional Linearized Inversion,” presented at 
the 57th Ann. Int. Mtg. and Expos., Soc. Explor. Geophys., New Orleans, 
Expanded Abstracts, 744-746. 

M. A. Hooshyar and M. Razavy (1983). ”A Method for Constructing Wave Ve¬ 
locity and Density Profiles from the Angular Dependence of the Reflection 
Coefficient,” J. Acoust. Soc. Am. 73, 19-23. 

M. A. Hooshyar and A. B. Weglein (1986). "Inversion of the Two-Dimensional 
SH Elastic Wave Equation for the Density and Shear Modulus,” J. Acoust. 
Soc. Am. 79, 1280-1283. 

P. M. Morse and H. Feshbach (1953). Methods of Theoretical Physics (McGraw- 
Hill, New York). 

S. J. Norton (1983). "Generation of Separate Density and Compressibility Images 
in Tissue,” Ultrasonic Imaging 5, 240-252. 

S. J. Norton (1987). "Three-Dimensional Seismic Inversion of Velocity- and 
Density-Dependent Reflectivity,” Geophys. J. R. Astr. Soc. 88, 393-415. 

A. Ozbek and B. C. Levy (1987). "Inversion of Parabolic and Paraboloidal Pro¬ 
jections,” Report LIDS-P-1665, Laboratory for Information and Decision Sys¬ 
tems, Mass. Inst, of Tech., Cambridge, MA; see also Proc. 21st Conf. Informa¬ 
tion Sciences and Systems, The Johns Hopkins Univ., Baltimore, MD, 36-42. 



Ozbek and Levy 


Velocity and Density Inversion 


30 


A. G. Ramm and A. B. Weglein (1984). ’’Inverse Scattering for Geophysical 
Problems. II. Inversion of Acoustical Data,” J. Math. Phys. 25 3231-3234. 

S. Raz (1981a). ’’Direct Reconstruction of Velocity and Density Profiles from 
Scattered Field Data,” Geophysics 46, 832-836. 

S. Raz (1981b). "Three-Dimensional Velocity Profile Inversion from Finite-Offset 
Scattering Data,” Geophysics 46, 837-842. 

V. T. Tatarski (1961). Wave Propagation in a Turbulent Medium (McGraw-Hill, 
New York). 

J. R. Taylor (1972). Scattering Theory (John Wiley, New York). 

A. E. Yagle and B. C. Levy (1984). "Application of the Schur Algorithm to the 
Inverse Problem for a Layered Acoustic Medium,” J. Acoust. Soc. Am. 76, 
301-308. 



Ozbek and Levy 


Velocity and Density Inversion 


31 


Figure Captions 

Fig. 1 2-D experimental geometry. 

Fig. 2a Scattering pattern due to velocity inhomogeneities. 

Fig. 2b Scattering pattern due to density inhomogeneities. 

Fig. 3 Coverage of Ur{J£) for a single array. 

Fig. 4 Frequency coverage of Ur{I^ and the "radiation pattern” of <Tnun(Af; fc; £i, £ 2 ) 
for the case when two probing waves are used, incident at right angles to each 
other. 

Fig. 5 Frequency coverage when of Ur{!^ and the "radiation pattern” of £ 1 , £ 2 ) 

for the case when two probing waves are used, incident at a 45° angle with 
respect to each other. 

Fig. 6 Velocity scattering potential model for the synthetic experiment. 

Fig. 7 Density scattering potential model for the synthetic experiment. 

Fig. 8 The backprejected image Ub"(x}, obtained assuming a constant density 
medium. 

Fig. 9 The inverted image obtained assuming a constant density medium. 

Fig. 10 The separate reconstruction of the velocity scattering potential. 

Fig. 11 The separate reconstruction of the density scattering potential. 
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